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FAST MODEL-BASED OPTICAL PROXIMITY CORRECTION 

TECHNICAL FIELD 

The present invention relates in general to manufacturing processes that require 
lithography and, in particular, to methods of designing photomasks and optimizing lithographic 
5 and etch processes used in microelectronics manufacturing. 

BACKGROUND INFORMATION 

Today's integrated circuits (ICs) contain features which are not easily resolved by 
available lithography tools, so proper printing of critical feature dimensions (CDs) requires that 

10 compensating adjustments be made to the IC shapes deployed on the mask. In so-called model- 
based optical-proximity-correction (MBOPC), the appropriate adjustments are determined by 
simulating the lithographic process, and in particular, providing a means to simulate the image at 
the wafer plane. Ultimately, a determination of the image formed in the resist layer (e.g. the 
latent resist image) is desired, but often the aerial image at the wafer plane is used as an 

15 approximation for the latent image in the resist. 

Conventional image simulation is typically done using the Hopkins integral for scalar 
partial coherent image formation, where the expression for the aerial image intensity / 0 is given 



where, 

h is the lens impulse response function (also known as the point spread function or PSF); 

j is the coherence; 

m is the mask transmission function; 



by, 



Equation 1 




25 



* indicates the complex conjugate; and 
r is the position of the image. 
The integration is typically performed over the mask. The expression: 
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Equation 2 

h(r -r')h*{jr- r') j(r' - r") 

is known as the Hopkins kernel, which is a fixed two-dimensional (2D) function for a given 
system. 

5 This 4-dimensional (4D) Hopkins integral (Equation 1) may be approximated as an 

incoherent sum of 2-dimensional (2D) coherent image integrals. This method of computing the 
Hopkins integral is known as the sum of coherent systems (SOCS) procedure. In the SOCS 
procedure, an optimal n-term approximation to the partially coherent Hopkins integral is: 

Equation 3 

10 I Q (r) = ±A k \(m®4> k Xrf , 

where ® represents the two-dimensional (2D) convolution operation, \,(/> k {f) represent the k th 

eigenvalue and eigenfunction of the Hopkins kernel, respectively, derived from the Mercer 
expansion of: 

Equation 4 

15 h{r'y{r'W-r') = £ Aj k (r% (?'), 

which suggests that a partially coherent imaging problem can be optimally approximated by a 
finite sum of coherent images obtained, for example, by linear convolution. Typically, the 
source and the mask polygons are decomposed (e.g. into grids or sectors), and each field image 
is computed as an incoherent sum of coherent sub-images (also referred to as component-images, 

20 or pre-images). The total intensity at an image point r in question is then the sum over all 

component images. In the SOCS approximation, the number of coherent sub-images that must 
be calculated is minimized, for example, by diagonalizing the image matrix to achieve an 
acceptable approximate matrix of minimal rank by eigenvalue decomposition. For example, 
even a large-fill source can be adequately approximated when the number of 2D convolutions n 

25 is about 10. 
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The mask transmission function, m , which can be approximated by a binary mask design 
pattern, usually consisting of polygons of one or more transmission, can be represented in 
different ways, such as grid cells or a sector decomposition, and the image at a point r may be 
represented by a finite incoherent sum of weighted coherent convolutions of the mask 
5 transmission function m and the eigenfunctions </> k . Each of the eigenfunctions, and its 

convolution with a possible mask sector, may be pre-computed, thus providing a fast method of 
computing the aerial image / 0 . The aerial image is often assumed to be an adequate 

approximation of the resist latent image. To account for resist effects, previous methods have 
used post processing the aerial image by applying a lumped parameter model or by applying 
10 resist blurring. However, in extending model-based-optical-proximity-correction (MBOPC) to 

the sub-lOOnm dimensions of next-generation IC products, the prior art has a number of 
limitations. 

For example, the scalar treatment is applicable for numerical apertures (NAs) less than 
about 0.7, where the angle between interfering orders is fairly small in resist, so that the electric 

15 fields in different beams are almost perpendicular or anti-perpendicular when they interfere. 

Under these circumstances, superposition is almost scalar. However, numerical apertures (NAs) 
of at least 0.85 must be employed when extending optical lithography to the sub- 100 nm 
dimensions. At the resulting steep obliquities the standard scalar Hopkins integral becomes 
inaccurate, and the vector character of the electric field must be considered. It is known that this 

20 can be accomplished by calculating independent images in each of the global Cartesian 

coordinates of the electric field (E x , E y , E z ), and then summing these images over each of the 
allowed orientations of the illuminating polarization, and again over each of the different source 
directions that illuminate the mask. However, the computational efficiency of this procedure is 
not adequate for MBOPC. 

25 Adam et al. (K. Adam, Y. Granik, A. Torres and N.B. Cobb, "Improved modeling 

performance with an adapted vectorial formulation of the Hopkins imaging equation," in SPIE 
vol. 5040, Optical Microlithography XVI, ed. Anthony Yen (2003), p. 78-91) and Hafeman et al. 
(in SPIE vol. 5040, Optical Microlithography XVI, ed. Anthony Yen (2003), p. 700) disclosed 
an extension of the scalar Hopkins imaging equation to include vectorial addition of the 
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electromagnetic (EM) field inside a film. However, the approaches of Adam et al. and Hafeman 
et al. ignore the z-component of the field, and do not take into account effects such as lens 
birefringence, tailored source polarization, or blur from the resist or mask. 

The scalar treatment of image formation ignores illumination polarization, and assumes 
5 that the lens and resist surfaces introduce negligible partial polarization. However, OPC for 

future lithography will also have to take into account the polarization properties of the lens itself. 
As NA and lens complexity increases, the polarization state of the output beam is changed by 
cumulative polarization-dependent reflection losses that arise at the surfaces of the lens elements. 
Attainable performance in the antireflection coatings that inhibit these losses will become 

10 increasingly poor as wavelength decreases into the deep ultraviolet (UV). Forthcoming 157nm 
lenses will be birefringent even in their bulk substrates, due to spatial dispersion in the element 
substrates. Skew incidence at beamsplitter and other coatings will distort the light polarization. 
There is also interest in sources that deliberately introduce polarization variations between rays 
in order to minimize the contrast loss that arises in transverse magnetic (TM) polarization at high 

15 NA (i.e. polarization-tailored sources). The exposed image is also impacted by the resist film 

stack, due to multiple polarization-dependent reflections between the various interfaces. The 
anti-reflective (AR) films that inhibit these reflections can be less effective over the broad 
angular ranges that arise at high NA. Refraction at the top surface of the resist gives rise to 
spherical aberration in the transmitted image. 

20 Finally, minimum feature sizes in next generation ICs are beginning to approach the 

resolution of the photoresist, and this must be taken into account during MBOPC. 

Methods are known for modeling each of these phenomena over mask areas of modest 
size. Vector images can be calculated by summing over all illumination and image-plane 
polarizations, and over all source points. 

25 Resist blur may occur due to the finite resolution of the resist. Propagation in spatially 

dispersive media has been treated in the physics literature, and this analysis has now been 
applied to 157nm lithographic lenses. Finite resist resolution can be treated in an approximate 
way by blurring the exposed optical image using a resist kernel, which in the frequency domain 
is equivalent to frequency filtering the 4D Hopkins integral. It is known that a post-exposure 
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blurring in the resist arises in the chemical image that is acted on by the resist developer; this 
resist blurring can be accounted for by convolution of the optical image with a blur function, or 
equivalently by attenuation of the image spatial frequency content by a modulation transfer 
function (for example, see J. Garofalo et al., "Reduction of ASIC Gate-Level line-end shortening 
5 by Mask Compensation/' in SPIE v.2440 - Optical/Laser Microlithographv VIIL ed. Timothy A. 

Brunner (SPIE, 1995), p. 171.). Hoffnagle et al. (in "Method of measuring the spatial resolution 
of a photoresist," Optics Letters 27, no.20 (2002): p. 1776.) have shown that the resist 
modulation transfer function (also known as the effective resist MTF) can be determined for a 
particular spatial frequency from critical dimension (CD) measurements in resist that has been 

10 exposed to a pair of interfering sine waves (this interference pattern exhibiting essentially 100% 
optical contrast). The latent image contrast deduced in this way can be substantially less than 
100%, and this modulation loss becomes steadily more important as feature-sizes come closer to 
the resolution limit of resists. For example, at a pitch of 225nm, the UV2 resist analyzed in 
Hoffnagle et al. transfers to the latent chemical image only about 50% of the modulation in the 

15 exposing optical image. More modest contrast losses arise even at relatively coarse spatial 

frequencies, and can give rise to proximity effects of relatively long range. 

However, such methods to account for resist blurring involve a direct spatial domain 
convolution of the resist blur function with the continuous optical image, which is relatively slow 
to compute, as compared to the polygons of the kind that are deployed on the mask, which are 

20 relatively simple to handle, and can be relatively fast to compute. Unfortunately, future MBOPC 
will only be practical if a way can be found to calculate these effects very quickly. 

Other aspects of the MBOPC process (e.g. etch simulation, simulation of mask 
electromagnetic effects, polygon edge adjustment, as well as many non-lithographic issues) are 
not directly considered in the present discussion; however, it is important to bear in mind that 

25 lithographic simulation is only one small aspect of the full IC design process. Mask polygons 

must be stored in a recognized data structure of complex hierarchical organization; the patterns 
are arrived at after considerable effort that involves a long succession of circuit and device 
design software. Additional processing software converts the post-OPC shapes to a format used 
in mask writing. Because of the complexity of this computer-aided design (CAD) process, it has 
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become customary to use suites of software design tools (including tools for MBOPC) to avoid 
problems when migrating the chip data from one stage of the design process to the next. 
Improved MBOPC methods must be compatible with this design flow. It is also desirable that 
improved methods for lithographic simulation be compatible with existing OPC programs. 
5 The SOCS method gives a fast and accurate enough algorithm for computing the scalar 

aerial image (AI). This would be enough to satisfy and even exceed any practical requirements 
if only the AI was needed as the final output of the algorithm. However, because of the presence 
of the resist development step, further computation is needed to get the resist image from the AI. 
Images obtained using the SOCS method are calculated by pre-storing the corresponding 

10 images of all possible semi-infinite sectors from which the mask polygons are composed. The 

images of polygonal mask features can therefore be calculated very rapidly. 

These pre-stored tables are based only on the scalar Hopkins integral. To take into 
account such phenomena as resist resolution, vector imaging, resist thin-film effects, illumination 
polarization, and lens birefringence, we need to obtain new tables that are able to reproduce the 

15 effects of these phenomena as if the phenomena arose from an incoherent sum of coherent 

images. 

Accordingly, there is a need for an efficient method of simulating images that more 
accurately includes non-scalar effects (i.e. "non-Hopkins" effects) including the vector electric 
field, polarization effects, and is applicable for computing a variety of images including an aerial 
20 image or a resist image. In addition, it would be desirable to implement such a method that can 

be incorporated into existing computer codes without significant restructuring of the code. 



25 
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SUMMARY OF THE INVENTION 

Bearing in mind the problems and deficiencies of the prior art, it is therefore object of the 
present invention to provide a fast method for computing images that takes into account non- 
scalar (i.e. "non-Hopkins") effects. 
5 It is a further objective of the present invention to provide a fast method of computing 

non-scalar images that can be incorporated into existing computer codes without significant 
restructuring of such codes. 

Still other objects and advantages of the invention will in part be obvious and will in part 
be apparent from the specification. 
10 In accordance with the present invention, an efficient method and system is provided for 

computing lithographic images that takes into account vector effects such as lens birefringence, 
resist stack effects and tailored source polarizations, and may also include blur effects of the 
mask and the resist. In accordance with the present invention, these effects are included by 
forming a generalized bilinear kernel, which can then be treated using a SOCS decomposition to 
15 allow rapid computation of an image that includes such non-scalar effects. A source intensity 

distribution s(k s ) is provided. The source can be characterized by an intensity distribution 
s{k s ), and may be characterized by independent polarizations m , each of which is characterized 

by a polarization map. A projection impulse response function h m is provided, which preferably 

includes the vector impulse response function of the lens including lens birefringence, but may 
20 also be generalized to include vector effects in the resist film stack. Furthermore, in the case of a 

source with tailored polarization, the projection impulse response function h may also include a 
source polarization map. A generalized bilinear kernel V{f' y r") is formed by forming a bilinear 

autocorrelation of the source intensity distribution s{k s ) with the vector impulse response 

function h{k',k s ). A generalized bilinear kernel can be formed to include mask blur. Resist 

25 blur may be included in a generalized bilinear kernel, for example, where averaging over the 

depth of the resist stack and/or a range of focus positions may be performed. Preferably, a SOCS 
decomposition of the generalized bilinear kernel is performed. The generalized bilinear kernel 
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V is formed so as to be independent of the mask transmission function m(r ). In another aspect 
of the present invention, the image l(r) at a point r is then computed by combining, or more 
particularly, performing a bilinear integration of the generalized bilinear kernel with the mask 
transmission function. The resulting image l(r) can then be used to perform model-based 
5 optical proximity correction (MBOPC). 

In a preferred embodiment of the present invention, the SOCS decomposition of the 
generalized bilinear kernel may be performed. First the image plane within the region of interest 
(ROI) integration domain is gridded. If possible, in accordance with the invention, the ROI 
integration domain is preferably folded according to the symmetry of the system. The 

10 generalized bilinear kernel (GBK) is computed and tabulated at the grid points of the image 

plane, in the folded ROI integration domain. In this preferred embodiment, the tabulated GBK 
values are remapped to a reduced basis, and then the eigenfunctions of the GBK are calculated in 
the reduced basis, and converted back to the original grid. If necessary, the dominant 
eigenfunctions may be iteratively refined against the tabulated GBK values, for example, as in 

1 5 the Lanczos method. Then, the convolutions of the dominant eigenfunctions with the possible 

mask polygon sectors are pre-computed. 

According to yet another aspect of the present invention, the computation of the image 
includes the step of decomposing the mask transmission function into an appropriate set of the 
possible polygon sectors. Then, for each of the dominant eigenfunctions, a pre-image is 

20 computed by forming the coherent sum of the contributions from the pre-computed convolutions 

of the dominant eigenfunctions with the mask polygon sectors within the ROI. The pre-image is 
weighted, preferably by the eigenvalue of the eigenfunction, and in particular, more preferably 
by the square root of the eigenvalue. Alternatively, the weights may be determined empirically. 
Finally, the image l(r) is formed from the incoherent sum of the weighted pre-images of all the 

25 dominant eigenfunctions. 

The present invention has the advantage that non-scalar effects such as vector effects 
including tailored source polarization, lens birefringence, and resist stack polarization, as well as 
blur in the mask or resist, can be incorporated efficiently in image calculation. 
The present invention provides the following advantages: 
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1. The avoidance of extremely large computational errors which can result from taking 
the derivatives of a function which already has been computed with rounding errors. 

2. Reduction of the time complexity of the computation of the resist image to that of the 
computation of just the aerial image. In accordance with the present invention, a generalized 

5 bilinear kernel is formed that has the same or similar form as the bilinear kernel used to form the 

aerial image according to the SOCS approximation of the Hopkins model. In other words the 
invention reduces the complexity of the resist image computation to the complexity of obtaining 
just the aerial image, or by at least a factor on the order of 10. 

3. Improved accuracy of the image due to the inclusion of non-scalar effects, including, 
10 but not limited to effects of the vector electric field, tailored source polarization, lens 

birefringence, defocus variations, resist stack and blur from the mask and the resist. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

The features of the invention believed to be novel and the elements characteristic of the 
invention are set forth with particularity in the appended claims. The figures are for illustration 
purposes only and are not drawn to scale. The invention itself, however, both as to organization 
5 and method of operation, may best be understood by reference to the detailed description which 

follows taken in conjunction with the accompanying drawings, in which: 
FIG. 1 illustrates a lithographic projection system. 

FIG. 2 illustrates a mask layout and a region of influence (ROI) according to the present 
invention. 

10 FIG. 3 illustrates the vector propagation of a ray through the projection system. 

FIGs. 4-6 illustrate flowcharts of the method according to the present invention. 

FIG. 7 illustrates a computer system configured to perform the method according to the present 

invention. 



10 



FIS9-2004-0011 



DETAILED DESCRIPTION 

In the following description, numerous specific details may be set forth to provide a 
thorough understanding of the present invention. However, it will be obvious to those skilled in 
the art that the present invention may be practiced without such specific details. In other 
5 instances, well-known features may have been shown in block diagram form in order not to 

obscure the present invention in unnecessary detail. 

Refer now to the drawings wherein depicted elements are not necessarily shown to scale 
and wherein like or similar elements are designated by the same reference numeral through the 
several views. 

10 FIG. 1 illustrates a schematic diagram of Kohler projection system 10 used in optical 

lithography. An energy source 100 provides actinic energy (actinic energy is that illumination 
energy that will drive the reaction in the resist), which may be projected through a condenser 
lens 1 10 to a patterned mask or reticle (120). Together, the source 100 and the condenser lens 
1 10 may be referred to as the illumination source 115. The light diffracted from the mask 120 

15 then passes through a projection lens 130 (having an entrance pupil 129 and an exit pupil 135), 
and is then projected onto a wafer 140, which typically includes a resist layer 142 and an 
underlying substrate 144. The resolution of the projection system 10 is related to the amount of 
light controlled by the lens aperture stop 125, as indicated by the maximum angle 6 of the 
projected light relative to the optical axis 150 along the z direction, illustrated in FIG. 1, where 

20 26 is the range of angles that converge from the lens 130 on to each image point. The 

numerical aperture of the projection system is given by NA = ns'm 6 , where n is the index of 
refraction of the medium surrounding the lens (which is about 1 for air). For the present 
discussion, the region 122 between the mask 120 and the projection lens 130 is hereinafter 
referred to as the object space, and the region 132 between the projection lens 130 and the wafer 

25 140 will be referred to hereinafter as the image space. 

Each image calculation from MBOPC programs typically involves only those mask 
shapes that are contained within a small region (the so-called Region Of Interest, or ROI) that is 
centered on the image point of interest F . For example, FIG. 2 illustrates a mask layout 300 
comprising mask polygons 301. The image point of interest r 141 has a corresponding point on 

11 
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the mask layout r\r" (indicated reference numeral 121). Only mask polygons 301 within a 
certain zone of influence from the image point is considered in the computation of the image 
intensity, which is referred to hereinafter as the "region of interest" or ROI 3 1 0. This ROI 310 
region is only modestly larger than the lens resolution, e.g. the ROI 310 might have a radius of 
5 l(Lim (at IX scale) for a lens used to print 90nm features. When a particular ROI field is 

evaluated, one does not calculate image intensity at points near the ROI boundary. Usually, only 
the intensity at the center point 121 of the ROI 310 is evaluated; alternatively, as disclosed in K. 
Lai et al., "Improvement of performance in sector-based OPC engine utilizing efficient polygon 
pinning method and simultaneous computation of multiple sample points," (co-assigned U.S. 

10 Patent Application Ser. No. 10/694,473, filed on October 27, 2003, the contents of which are 

hereby incorporated by reference in its entirety) the intensity may be evaluated at several interior 
points, so long as these points are far from the ROI boundary. Under this assumption, we can 
assume for purposes of discussion that periodic boundary conditions apply. However, this 
assumption need not be made during actual calculation of the sector tables. 

15 Suppose we then interpret the ROI 310 as the unit cell of a grating, where the size of the 

grating 300 is large compared to the ROI 310, but still microscopic on the scale of the projection 
lens 130. Referring to FIG. 3, thus, considering the coherent image component that arises during 
illumination from a single source point 1 1 1 (here the energy source 100 and condenser lens 1 10 
in FIG. 1 are considered together as the illumination source 1 15), the beams which diffract from 

20 the ROI 310 will have a subtended width (not shown) in object space 122 that is quite small on 

scales characteristic of the projection lens 130. Each such beam is focused into the entrance 
pupil 129 under Kohler illumination. As the diffracted order propagates from the entrance pupil 
129 to the exit pupil 135 and then to the image plane (e.g. at the wafer 140), the small width of 
the beams diffracted from the ROI 3 10 allows it to be treated it as a single ray being focused by 

25 the lens 130. 

Vector and Polarization Effects 
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Referring to FIG. 3, consider how the various lithographic effects described above can be 
accounted for during propagation of a diffracted order, represented by ray 201 from a point 121 
on the mask 120 through the projection lens 130 to a point r 141 on the image plane 142 (which 
is typically within the resist layer 142). A meridinal plane 220 can be defined as the plane 
5 containing the path of the ray 201 and the local optical axis 150' (which is the ray that passes 

through the center of the aperture stop 125, but passes locally through the mask point 121 and the 
focal point on the image plane r 141). Consider first the simple case in which the source 100 is 
either linear polarized or unpolarized, and where the projection lens 130 exhibits ideal, 
nonpolarizing behavior. Even with an ideal projection lens 130 and a source of fixed 

10 polarization and low illumination-NA, the polarization along a propagating ray 201 will undergo 
rotation of sorts, since transversality must be maintained as the ray changes direction within the 
projection lens 130. Since lithographic lenses are telecentric in both object and image spaces, we 
can consider the image point of interest F 141 to lie on the local optical axis 150" (taken to lie 
along z). Optionally, the effect of aberrations specific to a specific image location can be taken 

15 into account by applying the aberration in the aperture stop appropriate to the field location, as is 

known in the art. 

For a ray propagating in a direction the polarization of the 

ray can be tracked against an (s,p) basis, where the p component of the ray polarization lies 
within the meridional plane 220 and the s component of the ray polarization is orthogonal to the 
20 p component, where: 

Equation 5 

s ~T — T 

|zxk| in object space 122; 

/S /N 

p = sxtc 

Equation 6 
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s =-r: — 77- 

\z x k I in image space 132; and 

/? =s x/r 

Equation 7 

|z x *| in the resist layer 142, 
p -s x/r 

where the ray directions in image space 132 (£') and in the resist layer 142 (jf) are determined 
from the object space 122 direction (ft) using the magnification (sine condition), and the resist 
index (Snell's law). Note that we have the identity 
Equation 8 

s-s - s 

since the sign change between the s component of the ray polarization of Equation 5 and the s' 
component of Equation 6 simply represents the inversion of the ray direction by the lens 130 (i.e. 
the magnification is negative). Note also that a distinction is drawn between image space 132 
and the resist layer 142. Although it is not required by the present invention, preferably, the 
propagation from mask to wafer is treated using polarization ray tracing, but propagation into the 
resist stack is treated using a conventional plane wave thin-film approach. 

To make large-area lithographic simulations practical, it is customary to make 
simplifying assumptions about mask diffraction. This involves two approximations; first, the so- 
called thin mask approximation, and second, an assumption that the projection lens impulse 
response function (which arises from spatial invariance) is also invariant with respect to 
illumination direction. For the present invention it is necessary extend both approximations to 
the case of a polarized source. To some extent both approximations are improved at large 
reduction ratios, where propagation angles of collected orders are relatively small and mask 
features relatively large; however even at low mask-side NA both approximations can be 
invalidated by strong topographic effects in the patterned mask. 
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Vector Effects 

At low mask-side NA there is no difficulty in representing the source polarization 
5 distribution by a 2D map which specifies the (2D) polarization direction of each source point. 

However, the reduction to 2D becomes somewhat arbitrary when mask-side NA is non- 
negligible. By convention we will refer to the coordinate directions for the source map as x and 
y (i.e. we take the optical axis to lie along z ). The illuminating direction will be specified by 
the (global) x and y direction cosines, while local x and y coordinate directions will be used to 
10 specify the illuminating polarization along each ray, for example by defining the local x axis to 
be perpendicular to the plane 220 containing the ray and the global y axis (with the local y axis 
then being perpendicular to the x axis and the ray 201). The unit vectors along the x and y 
directions can be expressed as follows: 
Equation 9 

, = y-(y-K s )K s 

15 |y-(y-*i)*s| 

where x.y are the unit vectors of the global coordinate system, and k s is the illumination 
direction. 

Unpolarized or partially polarized sources are treated by independently superposing the 
system kernels calculated for independent polarized sources. 

20 The assumption of spatial invariance (i.e. the assumption that the projection lens 130 can 

be characterized be a fixed lens impulse response function within each local imaging field) 
primarily involves the projection lens 130. However, to treat partially coherent illumination an 
associated assumption is conventionally made about mask diffraction, namely that mask 
diffraction can properly be treated by applying the same fixed lens impulse response function to 

25 the mask shapes even as the illuminating direction shifts from source point to source point (with 

the amplitude transmitted through the mask shapes changing only by the illumination tilt phase). 

15 



FIS9-2004-0011 



The present treatment likewise requires that the projection lens exhibit a spatially invariant 
impulse response. It is not strictly necessary that this impulse response be invariant with respect 
to source direction, even for a single source polarization. However, the pre-calculation phase of 
the invention is significantly lengthened unless the impulse response obeys a constrained 
5 dependence on source direction, which is hereinafter referred to as a generalized source- 

invariance. The simplest such source-invariance model constitutes a natural extension of the 
assumption conventionally made in the scalar Hopkins treatment; the projection lens is found to 
be governed by two vector impulse response functions that are independent of illumination 
direction (except for a phase tilt), with each vector response corresponding to polarization along 

10 one of the two source polarization directions, x or y . Even though, at non-negligible mask- 

side NA, the polarization of the illuminating ray cannot really be completely invariant as ray 
direction is varied, since mask-side NA is reduced by a factor of at least 4, the polarization of the 
illuminating ray is almost invariant as ray direction is varied. 

Efficient pre-calculation is also possible with more complex models that do not require 

15 approximate source invariance. For example, the source 1 15 can be divided into independent 

regions (e.g. the poles of a quadrupole) that are each governed by a pair of source impulse 
response functions. In accordance with the present invention, efficient integration to determine a 
generalized bilinear kernel does not require full directional source invariance within the region; 
only that the response be a sum of separable functions of the illumination and collection 

20 directions (i.e. diffraction directions). 

The present invention uses the thin mask approximation, as known in the art, in that the 
amplitude of the diffracted light is calculated as the Fourier transform of the mask patterns (or, in 
some embodiments, as the Fourier transform of a blurred rendition of the mask patterns; in 
general the function must be linear). It has been shown that the thin-mask framework can be 

25 extended to accommodate mask topography effects by replacing the true mask patterns with 

effective patterns whose Fourier transform more closely reproduces the true diffraction spectrum. 
A simple approach along these lines is to bias each edge with a slightly increased chrome width. 

In a preferred embodiment of the present invention, consider the case of an ideal lens for 
a simple model of mask diffraction. Consider the on-axis source point to be polarized along a 
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direction E 0 that has no z component. The E 0 direction will usually have a component 

perpendicular to the meridional plane 220 of a given ray 201 ; we identify the projection of E 0 

perpendicular to the meridional plane 220 as the s component of the ray polarization, and 
assume initially that the amplitude of the s component is maintained in propagating through the 
5 lens 130 (this is true of an ideal lens; the general case will be considered later). In our initial 

simple scenario the s amplitude of the electric field is taken to be equal to this projection factor 
multiplied by the usual Fourier component of the mask diffraction pattern, which (due to 
relatively low NA) is calculated on a scalar basis (as is done in the prior art). The s component 
of the electric field is the given by, for an ideal lens: 

10 Equation 10 

5 V ' |zx£| 

where A{k) is the scalar order-amplitude, given in Equation 13 below. When the s component is 
mathematically subtracted from E 0 , the remainder is not transverse to the ray 220, and therefore 
does not represent a true p polarization component. Since we are treating diffraction from the 

15 mask 120 in an approximate way (given the relatively low-NA), we will initially assume that the 

magnitude of the p component can be calculated by finding the magnitude of this in-plane E 0 
component (multiplied by the appropriate scalar diffraction order amplitude A{k), The direction 
of this p component is of course taken to lie along the transverse direction p . We take the sign 
of the p component to be sign{p -E 0 ). 

20 After some algebra, we find that the p component is then given by: 

Equation 11 

E P p = [sign(p • E^E 0 - (E 0 - $)s|)a(£)p 
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where we have used z ■ E 0 = 0 . Note that Equation 1 1 will obey the identity 
Equation 12 

since the x and y components of k (the direction of propagation of the ray 220 in object space 
5 122) and K (the direction of propagation of the ray 220 in image space 132) only change in the 

ratio of the magnification, which cancels between numerator and denominator. 

In our model the electric field strength along the ray 220 will also be multiplied by the 
scalar diffraction efficiency A{k,k s ) of the mask 120 in the ray propagation direction k , which 

for unit illuminating intensity from source direction k s is given by: 
10 Equation 13 

A{k,k s )= Hdrm{r¥ { * s+ * } \ 

where m(r) is the mask transmittance function. 

In an ideal lens the s and p polarizations will not mix, and ray intensity will not be 
absorbed in the lens. As discussed- above, we can identify the ray 220 with a single order 

1 5 diffracted from the mask 1 20, truncated with a cross-section bounded by e.g. the ROI 3 1 0 in the 

mask plane (or some multiple thereof, taken here to be 1). The cross-section area subtended in 
the wafer plane is equal to ROI*M 2 (where M is the magnification of the lithographic system, 
and typically M = -0.25 ). However, this area is only subtended obliquely, i.e. it is 
foreshortened as seen from the ray direction, and the perpendicular cross-section area is therefore 

20 reduced by a factor y, where y denotes the z component of k , the direction of propagation in 

image space 132. This concentration causes an increase in the electric field strength. In an ideal 

lens the electric field along the ray (diffraction order) therefore changes by a ratio M yjy/y 

between object and image. We suppress the constant factor M for simplicity, and define the so- 
called obliquity factor as: 

25 Equation 14 
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We should also note that it is customary to describe mask-plane quantities in so-called "IX" 
wafer-scale coordinates, often without explicitly saying so. We will follow this custom below, 
referring to both the mask-scale and wafer-scale integration domains as "the ROI". 
5 The s component of the ray 220 is transverse to the resist plane of incidence (the resist 

plane of incidence is the plane containing the ray and the normal to the resist surface), so the s 
vector component propagates into the resist 142 with a transmittance T s (Az r ) , where T s (Az r ) is 

the coefficient for s transmission through the wafer film stack 142 to a depth Az r . Using 
standard thin film methods one can also calculate a transmission coefficient r p (Az r ) for the p 
10 component. The transmittance coefficients calculated by these standard methods provide the 
tangential component of the field (i.e. the component parallel to the interface of the resist stack 
142), so we define 

Equation 15 

1 5 as the transmission coefficient for the magnitude of the p amplitude. In the equations below we 

will suppress the Az dependence for brevity. 

If the bottom surface of the resist layer 142 is anti reflected, the electric field E" in the 
resist 142 will then be proportional to the quantity Q{k), that is: 

Equation 16 

20 E\k) = A{k; k s )0{k)Q"{k) , 

where Q{k) is defined as 
Equation 17 

Q(/c) = T s E s s"+f P E P p'. 
Note that the electric field in image space is similarly given by: 
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Equation 18 



E\K)=A{k'-k s )0^')Q'^') 



where the obliquity factor O^) in both the resist (Equation 16) and image space (Equation 18) 
is given by the same factor given by Equation 14. 

By making use of Equation 5, Equation 6, Equation 7, Equation 8, Equation 10, Equation 
1 1 , Equation 12 and Equation 15, we find after some algebra that inside the resist layer, with an 
index-matched substrate: 

Equation 19 



where k" is the projection of k" transverse to the local optical axis within the resist stack 142: 
Equation 20 



with the last form representing a vector of Cartesian components. Note that the first term t s E 0 
on the right of Equation 19 is essentially the standard, near-scalar result that obtains at low NA, 
while the remaining terms can be thought of as vector corrections. However, the r s factor 
appearing in the first term of Equation 19 can have a significant angle dependence at high NA. 
The term in square brackets Equation 19 is of order unity as k* — > 0 because r s f p , y* -> 1 at 

normal incidence; we therefore have Q^t s E q in the low-NA limit, as expected. 

The beam cone angle is more pronounced in air than in resist, so the image in air (the so- 
called aerial image) shows larger departures from scalar behavior at high NA than does the resist 
image. Thus, one must be cautious in approximating the resist image with the aerial image at 
NA's where vector interference effects are important. Bearing this caveat in mind, we note that 
Equation 19 provides as a special case a simple expression for the electric field E' along rays in 
the aerial image space 132: 

Equation 21 




k'=k"-{k'z)z = {a,p<d}, 



20 



FIS9-2004-0011 



e(*)=E 0 -%4W+rt]. 

Q"(k) of Equation 19 can be generalized to include the case of a stack of films 142 on 
the 144 wafer, e.g. antireflection layers above and below the resist, and/or other process films. 
Transfer into a general film stack 142 is not governed by a single transmittance T ; in addition to 
5 the transmitted down-traveling wave having amplitude u , the resist layer 142 will in general also 

contain an up-traveling wave having amplitude v that has reflected from the substrate 144. To 
solve for these amplitudes we divide the film stack 142 into an upper and lower substack, 
separated by the depth-plane within the resist at which we want to calculate the image. Standard 
thin-film methods can be used to calculate the up-traveling and down-traveling amplitudes, v,w , 
10 respectively; separate calculations must be made for the up and down-traveling s components 
(v 5 and u s , respectively) and the p components ( v p and u p , respectively). Thin-film methods 

conventionally provide transfer coefficients for the tangential field components, since the 
equations that are solved for these coefficients express the continuity of tangential components 
across interfaces. Our Equation 26 below recovers the full electric field from these components. 
15 To solve for the fields, we first calculate the standard thin-film characteristic matrix for 

each layer, defined as 

Equation 22 

( i \ 
cos B sin/? 

Y > 
^iYsxnfi cos/? j 

where J3 = 2md cos0/A , n is the film index of refraction, and d is the film thickness of the 
20 layer. Y is the layer admittance, defined as Y = ncos0 for s polarization, and Y = n/cos& for 

p polarization, with 9 the ray angle inside the film (calculated from Snell's law). The film 
matrices are multiplied together to form substack matrices M ] for the upper substack and M 2 
for the lower, and full stack matrix M Tot = M,M 2 . The following four linear equations are then 
solved for u and v (eliminating the supplementary unknowns p ToJ and T Tot ): 

25 Equation 23 
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[Y inc {l-p Tot )) 



= M 



To/ 



'Tot » 



= Af, 



t y 



'Tot ' 



W + V 

V 1 resist ^ 

repeating this procedure for s and p polarizations. The field inside the resist 142 is then 
obtained from: 



Equation 24 



Scalar term 



where u is defined as u = ^/y'* an d similarly for v. 

In the case of a non-ideal lens, the electric field in image space can be determined by 
polarization raytracing, as known in the art. The relationship between image-space field E' and 
object-space field E can be represented as multiplication by a Jones matrix J : 
10 Equation 25 

{E> s X p }=J$')-{E S9 E p } 

where the Jones matrix / may be used to include birefringence of the lens. Equation 25 can 
then be employed in the following generalized version of Equation 24: 

Equation 26 

15 £>(*')= (v 5 +w 5 )[£x£']e; +[v p (kY + zk")+u P ^ 

Equation 25 and Equation 26 assume that matrix / is calculated using the s, p basis of Equation 

5, Equation 6 and Equation 7. Thus, the electric field E' in image space obtained using Q(k') of 
Equation 26 includes the vector effects of the film stack. 

The s component of the electric field E s (of Equation 10) and p component of the 

20 electric field E p (of Equation 1 1 ) are written in terms of a single coherent illuminating 



22 



FIS9-2004-0O11 



polarization E 0 . Unpolarized illumination can be handled as an incoherent sum of coherent 
contributions from orthogonal initial choices of E 0 , as will be discussed further below. 

Customized Source Polarization 

5 

Another important generalization is that of a source having customized or tailored 
polarization. For example, as NA approaches or exceeds 1, the image quality of certain low-kl 
features can be enhanced by polarizing off-axis illumination in a tangential direction (i.e. in a 
direction perpendicular to the tilted plane of incidence). Formally, such customized polarization 

10 distributions can be handled by making E 0 (and therefore Q) functions of the projected 

illuminating direction k s , which is defined by analogy with Equation 19 as 
Equation 27 

k s =ft s A^s'z)z^{a s ,fi S9 0}. 

This also allows us to handle the polarization of light diffracted by the mask in a more general 
1 5 manner. 

Vector Image Based on a Generalized Bilinear Kernel 

In accordance with the present invention, we now obtain an expression for the overall 
20 image l(r) at a point F in the resist 142 that comprises a generalized bilinear kernel. To do so, 

we calculate the total electric field vector at each position r in the resist 142 (from each source 

point and from each independent illumination polarization), then calculate |e| , and then sum 

over all source points a and all independent illumination polarizations k s . (The sum over 
independent source polarizations, ranging from m = 1 to m = m M(lx , with m Max typically equal to 
25 1 or 2, allows us to account for unpolarized or partially polarized illumination.) We obtain: 

Equation 28 
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m=l 0 

1ft Max 

= Y,ljd 2 r'Hd 2 r"m(r'/R)m(7"/R) 

7 

jj d%-p(kfjo(kf}Q m (k-,k s 

^ J/J 2 r' JJd 2 r'm(r'//?)m* (r'/*) Jfrf a *,^fc ^ "£X(' " ^fc, )• £ (r - r";* s ) 



ROI 



As discussed further below, in accordance with the present invention, the image l(r ) of 
Equation 28 is now expressed in terms of a generalized bilinear kernel, V(r' 9 r") 9 given by: 

Equation 29 

which is dependent only on the source function s(£ 5 ) and the impulse response function of the 

lens h (which may also include information about the resist stack), but is independent of the 
mask transmittance function m(r). 

In the last line of Equation 28 we have written the mask patterns m(r) as if the reduction 
ratio R=l/Abs(M) were 1, following the common practice of describing masks in M 1X 
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dimensions". The integrals in Equation 28 should technically extend to infinity, but in practice 
they are restricted for computational reasons to a domain (such as ROI 310) modestly larger than 

the lens resolution, e.g. ± -4A/NA. p{k) represents the aberrated scalar pupil of the system, 

which may include defocus aberration. We have also introduced the symbol h to represent the 
5 electric field distribution in the image plane due to illuminating the projection, and is referred to 

as a lens vector impulse response, defined as 

Equation 30 

K (r ~ r% )m \jd 2 k'p(k')o(k% (*';*" 

NA 

where h m is the projection impulse response for a given source polarization m , which may also 

10 include lens birefringence, and may also include defocus aberration corresponding to a defocus 

position different from zero. h m does not include the effect of limited resist resolution, but in 

general includes the effect of multiple reflections within the resist stack. However, if we use the 

aerial image field Q(k) (as in Equation 21) instead of the more general expression for the image 

field of Equation 26, then the projection impulse response function h m of Equation 30 provides 

15 the aerial image impulse response of the lens without resist effects. Moreover, if we retain only 

the leading term in Equation 21, £ 0 , we recover the scalar aerial image impulse response in 

Equation 30 rather than a generalized impulse response that includes vector interference, resist 
film stack, lens birefringence, customized source polarization, and varying partial degree of 
polarization in the source. 

20 It should be noted that Q{k) of Equation 24 or Equation 26 only provide the field at a 

particular depth in the resist. This is consistent with the common practice of applying OPC to an 
image calculated within a single fixed image plane; however, in many cases it may be preferable 
to base OPC on the exposing image as averaged through the depth of the resist. This depth 
average can be approximated by an average over several planes, spaced by e.g. Az = 0. 1 5/1/ NA 2 . 

25 (One can average over lens focus in a similar way, for example to account for variation in 

positioning the wafer against the image beam; some of this variation will be truly cumulative in 
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modern scanned lithography tools, since it is incurred when the wafer is scanned through the 
imaging field.) If the total number of planes is W , then using w as an index to denote the 

particular Q value obtained for each resist plane (from Equation 24 or Equation 26), the bilinear 

terms appearing in the expression for the image intensity l(r) of Equation 28 are: 

Equation 31 

W NA NA vv=l 

which provide the image in resist averaged over the depth of the resist. Though Equation 3 1 
requires that the field be calculated in multiple planes, the actual use of depth-averaged images 
entails no additional computation once this pre-computation is completed (as long as the same 
total number of kernels is used). Note also that only a single sum over w should be made when 

h appears pairwise in bilinear terms; each copy of h or Q is indexed with the same value of w . 

For simplicity we denote the inner integrals in the expression for image intensity l(r) of 
Equation 28 with the symbol V , which is a 4D function in the x,y components of r' and f* , so 
that the image intensity l(r) of Equation 28, in accordance with the present invention, becomes 
a double convolution over the mask patterns: 
Equation 32 

/(F) = \\d 2 ?' \\d 2 ?"V{T', T')m{? - r')m{r - r") , 

ROI ROI 

where 
Equation 33 

V(r', r") . \\d%s{k s Y'^tk (r% )■ K (r';k s ), 

which, in accordance with the present invention, is a generalized bilinear kernel that is 
independent of the mask transmission function. Thus, in accordance with the present invention, 
the image intensity l(r) of Equation 32 can be viewed as a generalization of the scalar Hopkins 
integral, but which additionally includes the effects of vector diffraction, resist film stack, and 
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tailored source polarization. In addition, according to the present invention, this leads to a 
generalization of the SOCS method, which reduces the 4D Hopkins integral to a sum of squared 
2D convolutions. Additionally, if the polygonal shapes of IC mask patterns are exploited, for 
example, as described by K. Lai et al. (in co-assigned U.S. Patent Application 10/694,466, filed 
5 on October 27, 2003, the contents of which are hereby incorporated by reference in its entirety), 

then the convolution integrals over mask patterns can be calculated rapidly. Note that if the 
index of refraction of the resist is equal to 1, then l(r ) of Equation 32 is equivalent to the aerial 
image. 

10 Mask Blur 

However, the mask shapes m(r) will not be exactly polygonal, since mask-making tools 
have finite resolution, causing corner rounding in the pattern. Roughly speaking, we may 
approximate this loss in mask definition by convolution with a mask blur function b(r). If we 
15 denote the blurred mask patterns as m\ we have in this approximation the following replacement 

for the image intensity l(r) of Equation 32 that includes mask blur: 

Equation 34 

/(F) = \\d 2 7d 2 ry{rr)rn'{T - r')m* (r - F"), 

where 
20 Equation 35 

m'(r) = jjd 2 Rb(R)m(?-R). 

Unfortunately, the image intensity l(r ) including mask blur as in Equation 34 requires 
considerably more computation to evaluate than Equation 32, since the functions m are 
continuously varying functions, rather than binary polygons. 
25 However, we can generalize the intensity l(r) of Equation 28 to include the effect of the 

mask blur function b{r ) as follows: 
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10 



15 



Substituting the mask shapes m of Equation 35 into the intensity /(?) of Equation 34 

gives: 
Equation 36 

I{r) = ll\jd 2 rd 2 r"d 2 R'd 2 R"V{r\ r')b(R')b(R'}n(r - r'-R')m*(r - r"-R'). 
Making the change of variables 
Equation 37 

r^(r'-/?')l >S = r'+R', d 2 r'd 2 R' = d 2 Kd 2 s\ 
r = l(r'-R'\ Kr = r + R\ d 2 r"d 2 R" = d 2 r'd 2 s" 9 

gives 

Equation 38 

/(F) = Jjrf Wx^V V, K")m{r - *V(r - S), 

where 
Equation 39 



vV/)e \\d 2 s'd 2 rv 
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where V 



Y + J',— + J"' I is the generalized bilinear kernel without the mask (see Equation 33). 



It is convenient to modify Equation 38 and Equation 39 with additional changes of variables, to 
obtain 



20 



Equation 40 

with 

Equation 41 



I(r)= tfd 2 r'd 2 rV{r',r')m{r - r')m{? - r"), 



28 



FIS9-2004-0011 



V'{r\r') = tfd 2 K'd 2 K y{7\ F>(r'- K')>(r'-K") 

= jjd 2 rd>K"iid%s(ky'°^ 

Note that Equation 40 images the nominal polygonal mask patterns m(r), rather than the blurred 
mask m{f ) that is used in Equation 34. 

Thus, in accordance with the present invention, the image l(r) of Equation 40 includes 
5 the generalized bilinear kernel V'(r\r") of Equation 41 , which includes blurring from mask 

fabrication, and accounts for blurring imposed by the optics (and by multiple reflections within 
the film stack) using a generalized vector model (also accounting for the contrast loss at high-NA 
that can occur with interfering vector fields). 

However, resist blur may occur due to the finite resolution of the resist. As discussed 
10 above, this resist blurring can be accounted for by convolution of the optical image with a blur 

function, or equivalently by attenuation of the image spatial frequency content by a modulation 
transfer function. However, such methods, which involve direct convolution of the resist blur 
function with the image, are inefficient. 

15 Resist Blur 

In accordance with the present invention, the effects of resist blur can be included 
efficiently in the image, for example by determining an effective resist blur function from the 
measured resist modulation transfer function ( MTF ) by 2D Fourier transform (Hankel 
20 transform): 

Equation 42 

oo 

g w (?)= g 2D ^\)=2x\dvvJ 0 {27i\?\v)MTF{ V ), 

0 

where g 2 o(^) * s the exposure response of the resist, and more particularly, the exposure 
response at a plane in the resist stack structure, and J 0 is a Bessel function. The exposure 
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response of the resist is typically measured, or estimated, or provided by the resist manufacturer. 
The lithographically relevant image is determined from the optical image according to: 

Equation 43 

1^)= \\d 2 Rl{R)g 2D (r-R). 

ROI 

5 Unfortunately, the optical image /(/?) is continuously varying, making Equation 43 difficult to 

integrate in the rapid fashion required for OPC. However, in accordance with the present 
invention, we can account for resist resolution using a modified generalized bilinear imaging 
kernel V'that adds resist blur to the extended list of phenomena included in Equation 41 : 

Equation 44 

10 I u ,Jr)= \\d 2 7'd 2 r'V(7'rr')m{r-7'')m\7-r'), 

where the modified generalized bilinear kernel V"(r\r") in accordance with the present 

invention is 
Equation 45 

(7 

w=l m=\ ROI 

1 5 which includes resist blur averaged over depth of the resist stack and/or a range of focus 

positions. The generalized bilinear kernel V" includes a projection impulse response function h 
combined with the resist blur function g 2D . It also includes mask blur b and vector imaging 
effects. The combination of any one of the resist blur function, the mask blur function or other 
resist stack effects, is referred to hereinafter as an exposure response. If W is 1 , then Equation 

20 45 represents resist blur at a single resist plane. 

Note that the last term of the generalized bilinear kernel of Equation 45 given by 

Equation 46 

w=l m=l roi 
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represents the bilinear impulse response of the lithographic subsystem including the (scalar or 
vector) projection impulse response, mask blur and resist blur. Additionally, in accordance with 
the present invention, the subsystem impulse response can include the scalar or vector projection 
impulse response plus resist blur (without mask blur), or scalar or vector projection impulse 
5 response with mask blur (without resist blur). 

Even though the generalized bilinear kernel V' of Equation 45 would ordinarily make 

use of the generalized vector impulse response h m of the lens (Equation 30), we have noted 
above that the impulse response h m can instead be calculated using either the electric field term 

Q of Equation 21 or the leading term E Q of Equation 21 in order to obtain a generalized bilinear 

10 kernel involving either the aerial image or the scalar aerial image, respectively. 

Optionally, mask blur, resist blur, depth averaging and focus averaging may be added 

later. 

Transmission Cross-Correlation Coefficients (TCCs) 

15 

If the generalized bilinear kernel V*(r',r") of Equation 45 is Fourier-transformed, we 
obtain a generalization of what are known as the transmission cross-correlation coefficients 
(TCCs): 

Equation 47 

rn 

where each capital-letter variable denotes the Fourier transform of the spatial domain quantity 
that is represented by the same symbol in lower-case. Variable H is an exception to this 
convention; we use H to denote the circular lens pupil P (e.g. a circular tophat function, equal to 
25 the Fourier transform of the scalar projection impulse response function) combined with the 

obliquity factor 0, or more generally an aberrated scalar pupil combined with an obliquity factor 
and any defocus in the lens that may be present. 
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Application to Resist Blur Models 

The generalized bilinear kernel V" of Equation 45 includes the blurring in the optical 
5 exposure introduced by multiple reflections within the resist stack, the effect of finite resolution 

of the resist material itself, and (if averaged through depth) the effect of defocus within the resist 
layer, as well as the inherent resolution limits of a finite-NA lens imaging a vector field. 
However, it does not include dynamic effects that occur when the developer (which is the 
chemical used to develop the patterns of exposed resist) interacts with the exposed latent image 

10 to print the developed circuit feature, nor the similar effects that occur when the developed resist 

feature is transferred into an inorganic circuit film. These dynamic transfer effects are small, 
since the resist is strongly hard-limited, e.g. it is characterized by a power-law nonlinearity with 
the resist constrast, T, which is typically >10. Thus, resist development follows a constant 
threshold model to first order, and the goal that a given feature edge be printed at a specified 

15 target position is approximately met by requiring that the intensity at that position equal a 

reference intensity. This reference intensity is often chosen as the edge intensity of a particularly 
critical feature, since the development process is generally adjusted to print the most critical 
feature edge at nominal. 

However, in reality resist contrast T cannot be infinite, so the printed feature will usually 

20 be biased slightly away from the nominal image contour. This can be accounted for as a bias, or 
equivalently (for small biases) as an effective change in the edge intensity. In some development 
models, such as the Brunner-Ferguson model, this bias (or intensity change) is modeled as a 
function of the intensity and intensity slope at the nominal edge position, denoted f(l,dl/dx). 
(Without loss of generality we consider the right edge of a horizontally oriented feature, so that 

25 image slope can be represented as a positive derivative along x.) If / is expanded in a series, 
we have for the effective intensity at the edge: 
Equation 48 

GO = /(?)+ /(.,[/(?)- l(r ref )]+ / (2 )[/(r)-/(F^l 
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where / is shorthand for dl/dx , r ref refers to the edge of the reference feature, and where 

and / (2) represent small correction coefficients. 

Coefficients f ( i) and f (2) can be determined analytically, or by fitting to empirical data. 
The effective intensity can equivalently be expressed in terms of different parameters Co, Cj, and 

c 2 , as 

Equation 49 

h ff ( F ) = c 0 + l(r){\ + c x ) + c 2 /(r) 
If the generalized bilinear kernel V for intensity / is obtained by the methods described above 
according to the present invention, it can immediately be modified to incorporate Co and Ci by 
direct substitution. The kernel for / can be obtained in several ways, for example by 
differentiating an eigenvector expansion of V . We are then able to calculate the effective 
intensity l eff directly, rather than postcalculating it from separate calculations of intensity and 

intensity slope, thus speeding execution time. In accordance with the present invention, an 
effective generalized bilinear kernel is obtained as: 

Equation 50 

K # (^r)=c 0+ (i + c>(^^ 

where L represents the number of retained eigenfunctions *¥ of the generalized bilinear kernel 
V s \ 4* is the partial derivative of *F with respect to x . The effective generalized bilinear 
kernel V eff can then be diagonalized into a final set of eigenfunctions. 

Other empirical terms can likewise be added to V eff . Given an general functional 

dependence of AI on the behavior of / in the vicinity of an edge, we can expand this 
dependence in a functional Taylor series, and obtain a lowest order linear term that takes the 
form of an integration of / with a first-order kernel Q (i.e. the first order term of the functional 
Tayor series) 

Equation 51 
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A/(F) = c 3 \\d 2 r'Cl{r -r)l{r') 

ROI 

where we have included an explicit fitting constant C3 (which may also be absorbed into Q ). Q 
may be determined from fits to CD data; for example, CD data that contains many different 
assist feature combinations. Regularizing terms can be added during the fit to ensure smooth and 
monotonic behavior in Q . Once Q has been determined, the fit can be updated against later 
data by simply readjusting the linear coefficient C3. We then have 

Equation 52 

V eff (r ', r') = c 0 + (1 + c, r') + c 2 J>, (*,* (r')% (r ')+ 4>* (?')*, (r "))+ c 3 A V Q (r r ') 

where AV n is obtained by replacing the resist blur function g 2D in the generalized bilinear kernel 
V"(r' 9 r*) of Equation 45 with a corrected influence function g Q given by: 

Equation 53 

8a (r) = \\d'7W-r)g 2D {r). 

ROI 

(Experimental determinations of the influence function are likely to determine g Q directly.) 
Note that the (l + q ) and c 3 terms in Equation 52 can be combined as a total blur function 

Equation 54 

Seg ( ? ) = 0 + c i )g 20(^^8 cii?)- 

In many cases g Q will be a longer range function than g 2D . Formally, we can thus regard g L as 
either a long distance component of the resist blur function, or alternatively as a separate (often 
density-like) influence function. If readjustment in Ci and C3 is not contemplated, it may be 
desirable to carry out a one-time fitting in which (by iteration) the eigenfunctions of V eff itself 
are made to serve in the c 2 slope term of Equation 52. This is convenient if the tabulated kernels 
are to be input to pre-existing OPC software that does not support post-application of fitting 
kernels. On the other hand, if the long range character of g a forces an increased number of 
retained eigenvectors, one may prefer to leave the c 3 term separate from the main image 



34 



i 

FIS9-2004-0011 



generalized bilinear kernel, V" , adding it to the intensity only in cases where the most accurate 
OPC calculations are needed (e.g. avoiding it when checking for print-through or basic image 
polarity, or when carrying out initial iterations of OPC). 

5 Generalization of SOCS Decomposition 

In accordance with the present invention, the generalized bilinear kernel may be used to 
generalize the SOCS decomposition. The dominant eigenfunctions of the generalized bilinear 
kernel operator V"(or alternatively, the Fourier transform of the dominant eigenfunctions of 
10 operator T" which are then Fourier transformed), and their associated eigenvalues are calculated. 
The eigenfunctions and associated eigenvalues ju satisfy the equation 

Equation 55 

Jp 2 rV(r',r>(F') = /^(F')- 

ROI 

Fast Pre-Computation 

15 

Equation 55 may be solved for *F and ju using matrix eigendecomposition methods. 

Typically, only the largest 10 or 20 eigenvalues and associated eigenfunctions need be 
calculated. (We refer to these as the dominant eigenelements.) The eigenelements may be found 
by approximating the Equation 55 integral as a summation on a grid, with grid step of order 
20 02A/(NA(\ + G Max )) , with o Max denoting the tilt of the most obliquely incident illuminating 

direction used. We may also set a lower limit of, e.g. 0.4 on a Max , even in particular near- 
coherent exposure processes in which the maximum a value actually employed is somewhat 
lower. Therefore, the generalized bilinear kernel V may be approximated as a matrix, with 
variables F' and r" stepped out along rows and columns respectively. (Though each of these 
25 variables is actually two-dimensional, their associated 2D grids are unwrapped into an unraveled 

linear grid along a matrix axis.) 
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Variable Grid Pre-Computation 

Matrix size can be reduced without sacrificing accuracy by using a variable-grid 
algorithm, in which a coarser grid is used where the generalized bilinear kernel is small and/or 
slowly varying. Let the x and y grid sizes at the kth row or column be Ax k and Ay k . If we 

then calculate the eigenvectors .,y .) of the modified discrete matrix equation 
Equation 56 

tw'( x j >yjh Z^fo » y* • x i ■ yj M*.- > ^ ) ■ 

where the regridded matrix V is given by 
Equation 57 

V(x s , y f , Xj ,yj)= V(x i , y f . , Xj , y y ^ Ax,- Ay,. Ax . Ay . , 
then the desired eigenvectors *F(x.,y .) of the imaging kernel V will be given by 

Equation 58 

^/ Ax-Ay. 

Reduced Basis Pre-Computation 

The dominant eigenvectors and eigenvalues (Equation 58) of the kernel matrix (Equation 
57) can be obtained by standard methods, such as the Lanzos method, which can avoid 
calculation of most non-dominant eigenelements, improving efficiency. 

Matrix size can be significantly reduced by remapping to a suitable reduced basis. 
Preferably, a suitable reduced basis may be obtained so that the eigenvectors of the original 
kernel matrix are approximated by the new reduced basis vectors, for example, reduced basis 
vectors from previously-calculated eigenvectors obtained under similar system parameters, or 
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reduced basis vectors formed from a coarse-grid solution. Similarly, eigenanalysis can be 
performed on such a reduced-basis matrix in order to obtain a good starting solution for the 
iterative eigensolution algorithm. To do so, the eigenvectors of the reduced-basis matrix are 
converted back to the original grid basis. Then, the eigenvectors for the original kernel matrix 
may be iteratively refined, if necessary, by methods such as the Lanzos method. 

The total element count in the kernel matrix (whether V , the generalized bilinear kernel 
in the spatial domain, or T , which is the generalized bilinear kernel in the Fourier domain) 
scales as the 4th power of the ROI, and the time required to calculate the eigenelements usually 
increases even more rapidly (for example, eigenanalysis computation time may scale as the 3/2 
power of the element count). In cases where the ROI must be unusually large (for example, 
when flare or other interaction from distant features must be included, or when the effective 
resist blur function has a long tail), it is particularly desirable to minimize the precomputation 
involved in calculating the kernel and its eigenfunctions. This is also important when many 
different system models must be computed, e.g. in data fitting. 

Application of Symmetry to Pre-Computation 

Evaluation time of the generalized bilinear kernel, in accordance with the present 
invention, can be reduced by a factor of 2 if the Hermitian symmetry of the kernel is exploited, 
i.e. if one exploits 

Equation 59 

V"{r\r") = V"*{r\r') 
to eliminate half the kernel evaluations. 

In many cases it is acceptable to regard the ROI boundary as circular; if so, both kernel 
evaluation time and matrix size can be reduced by a factor (7i/4) 2 by using a circular rather than a 
square domain (in each variable of the kernel). 

Both kernel evaluation time and eigenanalysis time can usually be reduced substantially 
by exploiting system symmetry. For example, projection lenses have nominal rotational 
symmetry if residual imperfections are neglected (at least, the underlying design form has a 
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nominal axis of rotational symmetry even when the physical lithographic lens itself uses an off- 
axis imaging field), and usually it is not desirable to include residual asymmetries when carrying 
out OPC correction, since aberration-specific correction would entail the use of tool-specific 
masks, and would also force single circuit elements to undergo multiple OPC solutions when 
5 repeated in different places across the chip. Defocus is in some cases an exception to this rule, 

i.e. defocus can be the one aberration which should be taken into account during OPC; however, 
defocus preserves rotational symmetry. 

Even though the present invention can handle systems that are not symmetric (and can 
treat images which are averaged over a scan path by using a single averaged 4D kernel), 

10 precomputation is much faster when symmetry is exploited. 

System symmetry also requires symmetry in the source. Residual source asymmetries, 
like residual lens asymmetries, are most often deliberately neglected during OPC; however, it is 
often the case that even the nominal source will inherently exhibit a lesser degree of symmetry 
than the lens. (For example, the source may have a dipole shape even when the lens is 

15 rotationally symmetric.) In virtually all cases, the (nominal) source shape will at least have 

bilateral symmetry about the x and y axes. In most cases the source will also be symmetric about 
the 45°, 135° diagonals (since this provides matching performance in vertically and horizontally 
oriented features), the principal exception being dipole sources. In some cases the source shape 
will have rotational symmetry as well (e.g. disk or annular sources). These shape symmetries are 

20 preserved in the source as a whole when the source is unpolarized, and should preferably be 

maintained when source polarization is deliberately customized to a more complicated 
distribution. Any of these customary source symmetries will be shared by the nominal 
projection lens (even with defocus), as well as by the isotropic films in the resist stack, and also 
by the isotropic mask and resist blur functions; thus it is the source symmetry which usually 

25 determines the symmetry of the entire system. 

Application of Bilateral Symmetry in Pre-Computation 
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In nearly all cases of interest for OPC the system will exhibit bilateral symmetry about 
the x and y axes. We will refer to this as "dipole symmetry", since this symmetry is obeyed 
even by dipole sources (which are perhaps the least symmetric nominal sources used in 
lithography). Note that a source with dipole symmetry actually has a four-fold mirror symmetry, 
5 i.e. 

Equation 60 

V(x 9 y ;x ,y ) = V[-x 9 y \-x ,y ), 
V(x\ y'\x", y") s V(x'-y';x"-y") . 
Equation 60 dictates a similar symmetry in the eigenfunctions of the generalized bilinear kernel 
10 V . To see this, note that if x ¥(x\y / ) is such an eigenfunction, then application of Equation 60 

after direct substitution of e.g. *¥(-x\y') into Equation 55 demonstrates that ¥(-*',/) will 
also be an eigenfunction of V with the same eigenvalue. We thus see by forming the linear 
combinations [*F (jc', y')± *F(- x\ y')]/V2 that the eigenfunctions of V will be either odd or even 
in x (and likewise y ), or, in the case of two-fold degeneracy in ju , that the associated 
15 eigenfunctions will be spanned by subspace eigenbasis functions that are odd or even only. 

Let us then specify the odd or even symmetry of along x with the parameter % x , so 

that g x = -1 when *F has odd symmetry, and +1 with even symmetry. Similarly, we use 
parameter % y to specify the y -symmetry of *F . By making the appropriate changes of variables 
in the Equation 55 integral, we can then remap the domain of integration to the positive x\y 
20 quadrant: 
Equation 61 



0 0 L 



^v{ x ',-y;xy) + ^y{-x\-y-yy) 



The integration domain in Equation 61 has 1/4 the area of that in Equation 55, and it need only 
be evaluated over 1/4 as large a range in x*,y* (namely 0< x" < ROI 9 0< y" < ROI ), since *F 
25 is determined by symmetry in the other (negative) quadrants (according to £ x >€ y ). The matrix 
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obtained by discretizing Equation 61 thus contains 1/16 the number of elements as would be 
obtained by simply discretizing Equation 55, and the eigenelements of an Equation 61 matrix can 
be found far more rapidly (e.g. in 1/64 the time if eigensolution speed scales as the 3/2 power of 
element count). However, Equation 61 must be solved four times (since there are four different 
combinations of ^ x ,^ y ), so in a typical case the total eigensolution time may improve by about 

16X when bilateral system symmetry is exploited in this way. Speed also improves because 
generalized bilinear kernel V of the present invention need only be evaluated at 1/4 as many 
points. Once the separate (dominant) eigensolutions are obtained for each combination of £ 

and £ v , they are merged and resorted to obtain the dominant eigenelements of V . Note that 

different ^,£ v combination generally produce different eigenvalues. 

Application of Quadrupole Symmetry to Pre-Computation 

Besides having symmetry about the x and y axes, the source (and the projection system 
as a whole) will in most cases also have bilateral symmetry about the ±45° diagonals. We will 
refer to this as "quadrupole symmetry", since it is obeyed by standard quadrupole source shapes; 
note that it implies an eight-fold mirror symmetry of the system about the boundaries of each 
octant. The generalized bilinear kernel for such systems will obey the relation 

Equation 62 

V(x\ y\x\ /) = V(y\x\ y\ x") 
while continuing to obey the Equation 59 and Equation 60 dipole and Hermitian symmetries. 
We can exploit the Equation 62 symmetry in somewhat the same way as was done above with 
dipole symmetry, but the quadrupole case is more complicated. If with a quad-symmetric kernel 
we substitute the transpose y ¥{y\x'') of an eigenfunction satisfying Equation 61 into the 
Equation 61 operator (i.e. we apply the folded dipole-symmetry kernel to the transposed 
[diagonally flipped] eigenfunction), then make a change of variables (swap) x — > y\y -> x in 

the integral, and then apply Equation 62, we see that (with appropriate choice of £ x and £ v ) 
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*F(y',x') will be an eigenfunction of V with the same eigenvalue as *P(x', /). Moreover, if 
£ =£ y , the transposed eigenfunction will have the same dipole symmetry as *F(x',;y'), i.e. it 
will satisfy Equation 61 under the same values of g x and £ . Thus, we first consider 
determination of the subset of quad eigenfunctions which obey the dipole symmetry £ = +1 or 
- 1 , where £ is defined by ^ = ^ x = £ y . This set of quad eigenfunctions is obtained by reducing 
Equation 61 to an integration over one octant: 
Equation 63 

ROl ROl 

Mr") = J dx \dy^{x\ y'Mx', /; ?')+ V(- x,-y'; ?")+ {{v{- x, /; ?*) + V(x',-y'; ?")} 

+ £45 (y(y > x '> r")+V{- y'-x'\ r") + £{v(- >', x'; r") + V(y',-x'; r ")})] 

where £ 45 can be +1 or-1. If Equation 63 is solved by reduction to a discrete grid (i.e. to matrix 
equations), its solution for given £ entails eigenanalysis of two matrix choices (corresponding to 
the choices £ 45 = + 1 and £ 45 = -1 ), each reduced in size 4X from the dipole symmetry case. If 

solution time scales as the 3/2 power of element count, this provides a net 4X speed 
improvement over the corresponding portion of the dipole eigenanalysis. 2X fewer kernel 
evaluations are needed. 

The remaining quad-symmetry eigenfunctions are those associated with g x = (i.e. 

they are solutions with what may be termed opposed or x- y -asymmetric dipole symmetry); 
these do not obey an octagonal folding of the kind exploited in Equation 63. This may be seen 
by substituting *¥(y\x') into Equation 61 as above; we find that in this case the transposed *F is 
an eigenfunction with the same eigenvalue as y'), but with opposite dipole parameters (i.e. 
with £ x -> £ v , £ v — » £ . Thus, to exploit quadrupole symmetry in calculating such 

eigenelements, we first solve a single opposed-symmetry combination using the dipole- 
symmetry folded kernel of Equation 61 (e.g. we choose g x = + l,£ v = -1 , and solve that case as a 

dipole). After doing so, the eigenfunctions for the remaining case (e.g. £ v = +!,£=-!) are 
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recovered immediately by taking the transpose of the first set (e.g. ¥(/,*')). This represents a 
net 2X speed improvement over the corresponding steps of the dipole solution. Likewise, only 
half as many kernel evaluations are needed if Equation 62 is exploited. 

Though quadrupole symmetry implies that ^(y^x') will be an eigenfunction whenever 
5 ^(jc', /) is an eigenfunction, this is a rather trivial result when % x = g = ±1 , because in those 

cases the octagonal folding symmetry of Equation 63 indicates that *F(y',jt') and *¥(x',y') will 
actually be the same function (to within multiplication by -1). However, the eigenfunctions for 
the remaining two dipole cases ( £ x = -£ v = +1 and £ = -g = -1 ), though related by a simple 

transpose, are distinct functions (with the same eigenvalue); thus the (equal) eigenvalues 
10 obtained in these cases are two-fold degenerate. When £ x = -g y , any linear combination of 

and ^(x'^y') will also be an eigenfunction having the same eigenvalue. 
Application of Rotational Symmetry in Pre-Computation 

15 If the optical system is rotationally symmetric, we can find the eigenelements using a 

numerical version of the method disclosed in R.M. von Bunau, Y.C. Pati, and Y.-T. Wang, 
"Optimal coherent decompositions for radially symmmetric optical systems," J. Vac. Sci. 
Technol. B 15, no. 6 (1997), p. 2412, hereinafter referred to as Von Bunau, the contents of which 
are hereby incorporated by reference in its entirety. (Von Bunau also found analytic expressions 

20 for the eigenfunctions of the focused Hopkins kernel in the case of disk or annular sources.) In 

all likelihood the majority of lithographic sources in use today are disk-shaped or annular, and 
these sources will be circularly symmetric if unpolarized or tangentially polarized (or if treated 
as scalar). Rotational symmetry is preserved under defocus, and also in the presence of isotropic 
mask or resist blur, and likewise within isotropic layers of the resist film stack. Rotational 

25 symmetry is a continuous symmetry, and as shown below this allows us to reduce the 

dimensionality of the eigenprobiem by 1, providing a considerable time saving (increasingly so 
as problem size increases). Moreover, since we are only interested in the dominant 
eigenelements, we can reduce the dimensionality from 4 to "2.5", in the sense that decomposition 
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of several 2D problems can replace decomposition of the symmetry-derived 3D kernel; in cases 
of practical interest the additional time savings is appreciable. Note also that with radial 
symmetry, all Fourier transforms are efficiently carried out using fast Hankel transforms. 
To see how these savings are realized, we consider the generalized bilinear kernel 
5 V'(r' 9 r") of Equation 45 in the case of an unpolarized source. In that case m Max is 2, and the 

sum extends over uncorrected images under x-polarized and y-polarized illumination; thus 
under the ideal lens model for the s polarization component of the electric field of Equation 10 
and p component of the electric field of Equation 1 1, in accordance with the present invention, 
the generalized bilinear kernel V" becomes: 
10 Equation 64 

ROI 

where j(r) is the coherence function occurring in the Hopkins kernel of Equation 2, and where 
the projection impulse response function is 
Equation 65 

15 h m (r)^\\d 2 Kh m (K)b(r-K). 

We therefore form over one octant of F the function 
Equation 66 

f{r,s)= j(2s)\hXr~sx+R}h^ 

and then fill in the remainder of the r domain by symmetry. After convolving / with g 2D by 
20 FFT, then, in accordance with the present invention, the symmetric generalized bilinear kernel 
V 0 can be expressed as a 3D function in polar coordinates: 

Equation 67 

V'(r',F')=VV/;r , /) = V ; '(r' ) r'/-^= f(r'cos{0)+r"ys\n{0),jr' 2 +r" 2 -2rV'cos(0)) 
where 6 = 6' -9" . 

25 If we then denote the m th order azimuthal Fourier component of V" as v m , i.e. 
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Equation 68 

we can reduce the eigenanalysis problem for a given value of m to a purely radial one (ID in 
each variable). When calculating those eigenfunctions of V" that are associated with the m th 
azimuthal order, we find that the j th such eigenfunction is given by 

Equation 69 

M = -7= ¥ ijm {r)cosm0 if m > 0, 

= -7=^.o( r )if m = 0, 

= -p-^y_ m (r)sinm0if m < 0 
Vr 

where ^. m denotes the j th eigenfunction of the symmetrized 2D radial kernel obtained from 

m 7 

Equation 70 

ROl 

Aj, m <P'jJr")= \dr'J7? v m (r\r")<p' jm (r'). 

o 

The eigenvalue associated with x V jm is X j m . Note that the eigenvalues are degenerate when 
m * 0 ; the eigenfunctions y ¥ j m defined in Equation 69 for + m and - m then provide a pair that 
span the subspace. A reasonable value for m is about 10. 

Application of the Generalized Bilinear Kernel to OPC 

Once the eigenelements are calculated, the image intensity l(r ) in accordance with the 
present invention may be calculated at OPC fragmentation points in much the same way as is 
done in conventional OPC methods. For example, this may be done by pre-storing a set of sector 
convolutions for each of the dominant eigenfunctions *F (for Manhattan geometries, this set 
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need contain only a single such sector convolution, representing convolution with a 90° corner; 
however other geometries require additional convolutions for other corner angles). The image 
intensity at point r is then calculated by a sum over M dominant eigenelements and L corners 
enclosed within the ROI: 



where r ( is the position of the / th corner, and where 8 { is +1 or -1 , depending on the order in 
which the corner occurs when tracing the perimeter of the polygon that contains it (see K. Lai et 
al., co-assigned U.S. Patent Application 10/694,466). The C Y tables can then be used by 
MBOPC software (including software that is designed for use with the Hopkins model). 

Equation 71 is easy to differentiate with respect to position r , for example in a direction 
perpendicular to the polygon's edge, and one can prestore tables of the derivatives of the . 
This means, as discussed further below, that resist development models which involve image 
slope can also exploit the fast computation speed provided by the present invention. 

At this point it is appropriate to consider some computational aspects of Equation 45 for 
the generalized bilinear kernel V'and Equation 47 for the transmission cross-correlation 

coefficient T* . Their evaluation is complicated by the dependence of h m on both r' and k s . 

(Of course, this evaluation is only required during the pre-calculation steps of the invention; once 
the eigenfunctions of the resulting kernel are obtained, Equation 45 imposes no computational 
burden when the invention is actually applied to the mask shapes of interest.) As noted in the 

discussion surrounding Equation 27 regarding a polarized source, the dependence of h m on k s 

arises with sources that have customized (non-uniform) polarization. The most well-known 
examples are probably annular sources in which the polarization is made tangential in order to 
image critical spatial frequencies (at high NA) in transverse electric field (TE) polarization. The 

polarization changes in the ray polarization E that are incurred between mask and wafer will 

depend only on the diffracted direction k" , but the polarization that diffracts from the mask will 



Equation 71 



M 



L 
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also be dependent on the polarization of the illuminating ray, and hence on k s . Let x x anc * X y 
denote the x and y polarization components of the illuminating ray (in the local coordinate 
system for source polarization described above), so that % x and x y provide a map of the 

polarization distribution of the source, which may be tailored or customized. Then, using the 
thin-mask approximation, the following expression for the p component of the ray polarization 
after diffraction from the mask is obtained: 
Equation 72 

E p (t,k s ) = A{k + k S )S (k s Jz, (*>,, y (ic,k ,)+ Xy {k s )D p y (r,k s )]p . 

(A similar expression holds for the s component.) Here A(fc, k s ) is a Fourier component of the 

mask pattern, given by Equation 13. We have chosen to break out the intensity s(k s ) along the 

illuminating ray as a separate term, making % x X y * e components of a unit vector. The D 

coefficients represent the dependence of diffracted polarization on illuminating polarization. For 
example, under the generalized assumption of spatial invariance discussed above, we would take 
the diffracted polarization in the case of e.g. an x' -polarized source to be independent of 

illuminating ray direction. We indicate this k s -independence by adding a superscipt ( H ) (to 

denote pupil H rather than source S ). If we further apply the specific model used above to 
derive Equation 10 and Equation 1 1, we have (see Equation 11) 

Equation 73 




(We have also absorbed the p unit vector into D^ H \ making it a vector quantity.) 

Rapid evaluation of the integrals comprising the generalized bilinear kernel V" of 
Equation 45 does not force the use of the simple spatial invariance model of Equation 73. More 
generally, we need only assume that D p x > (along with D p y 9 D s / , and D s y > ) can be written in 
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separable form within each region of the source. For example if we divide the source into a total 
of J regions, each referenced by index j , then we require that 

Equation 74 

for each of the J source regions. (Either D (H) or Z) (5) may be a vector quantity.) Noting that 
both the lens parameters in Equation 25 and the film-stack parameters in Equation 26 are 
independent of illumination direction (and bearing in mind the simple geometrical 
correspondence between the various ray direction variables k", k' and k , we can collapse all 

dependence in Q on ray direction into generalized D coefficients [which we denote 

that are specific to the source polarization (i.e. to the specific x and y components along a 

given illuminating ray) and to the source region j . The projection impulse response function 

h m Equation 30 then becomes 

Equation 75 

NA LV J 

■ Km ir - r')S py (k s )+ % m ^.{r - r')S„ y (k s ) 
+ Kj,,Ar - r% y fc )+ l n ,. s .Ar ~ r')S s , y \k s ) 

where we have introduced impulse response components h that do not depend on the 
illuminating ray direction (except through index j ), and where we have combined the 

illumination-dependent terms D (5) and S into single factors denoted S . 
The generalized bilinear kernel of Equation 45 becomes 

Equation 76 
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H a 

x 2[X ,y. P y fc'fcy fe )+ K ,. p/ (k')s Pi/ )+ ^(K'^fcJ+^.^.^K'^fc)] 

x J|rf 2 ^ 2D ^|)b(F' + ^-K'>(F' + ^-K') 

If the dot products are expanded out, the integrals for the generalized bilinear kernel V*(r\r") 

reduce to a sum of bilinear integral terms in h . Several approaches are available to evaluate 
these bilinear terms. For example, one method first considers the preliminary case of no blurring 
from resist or mask (Equation 33). In the no-blur case, V is given by a sum of terms having the 
form 

Equation 77 

v e , f;no . blur = l(n l;{r)\\d%s{k s %{k s )s; 

a 

=lXryl;{r)j eJ {r'-?') 

where subscripts e and / can refer to any terms in the sum, and where a coherence term J has 
been introduced to denote the result of the source integral, and where in the last line we have 
approximated the unblurred V term by a bilinear sum of its dominant eigenfunctions (which are 
denoted ©). 

When resist and/or mask blur is added, the expression for the generalized bilinear kernel 
V 0 (Equation 76, where the double prime on V denotes the presence of blur) contains blurred 

versions of the same set of bilinear terms in h that arise without blur (given in Equation 77). 
When blur is added, these terms take the form 

Equation 78 
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Substituting from Equation 77, this becomes 
Equation 79 

K f =2X,/ J^ 2 ^ 2D ^|k,, / (r , +^V;.e,/(?'+^) 

n=l 

5 where 

Equation 80 

w n^,f(r)= l\d 2 K.e n<ef (K}?(r-K). 

The functions in the generalized bilinear kernel V" of Equation 79 are smoothly varying and 
hence easy to integrate numerically in the spatial domain. The number of integrations involved 

10 is similar to that involved in the calculation of scalar TCCs in an ordinary Hopkins model 

(though of course many such e y f combinations must be evaluated in the most general case). 

Another approach to calculating V 0 is to use integration in the frequency domain {T* of 
Equation 47), which can be easier than spatial-domain integration (except in the case of simple 
convolution) if one explicitly determines the boundary intersections of the pupil H and source 

15 S functions. The spatial domain kernels V" can then be determined by inverse transform. 

Efficient integration bounds in the frequency domain are defined by the bandlimited system 
response, and the periodic repeats of the object structure that are implicitly generated by 
integration against a discrete grid in the frequency domain may mimic excluded mask structure 
(i.e. exterior to the ROI) in a more realistic way than does simple spatial domain truncation. 

20 More preferably, another approach is to employ a mixed strategy. If the Equation 77 

blur-free eigenfunctions © n e f are Fourier transformed (we denote their transforms by @ n e f ), 
the transform of Equation 79 becomes 

Equation 81 

^=1 
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which after tabulation can be inverse transformed to obtain V" f . The eigenfunctions may also be 

transformed. The advantage of this approach is that, as shown by the middle line of Equation 77, 
tabulation of the unblurred kernel V does not require repeated integral evaluations. (One-time 
Fourier transforms are required, but these precede the tabulation loops, and may be carried out 
5 quickly by Fast Fourier Transform (FFT).) Once the unblurred kernel is calculated and 

eigendecomposed, it can be converted rapidly to the frequency domain by FFT, where 
application of blur is again a simple matter of function evaluation (no integration). However, 
this procedure does require two eigenanalysis steps. 

In summary, the present invention is directed to an efficient method and system for 

10 computing lithographic images that takes into account vector effects such as lens birefringence, 
resist stack effects and tailored source polarizations, and may also include blur effects of the 
mask and the resist. In accordance with the present invention, these effects are included by 
forming a generalized bilinear kernel, which can then be treated using a SOCS decomposition to 
allow rapid computation of an image that includes such non-scalar effects. Referring to FIG. 4, a 

15 preferred embodiment of the method according to the invention includes providing a source 

intensity distribution s(k s ) (Block 510). The source can be characterized by an intensity 

distribution s[k s ), and may be characterized by independent polarizations m , each of which is 
characterized by a polarization map (such as X X *X\ m reference to Equation 72). A projection 

impulse response function h m is provided (Block 520), which preferably includes the vector 

20 impulse response function of the lens including lens birefringence (e.g. Equation 30), but may 

also be generalized to include vector effects in the resist film stack as in Equation 31. 
Furthermore, in the case of a source with tailored polarization, the projection impulse response 

function h may also include the source polarization map (% x 9 % ), for example, as illustrated by 
Equation 75. A generalized bilinear kernel V(r\r") is formed (Block 530) by forming a bilinear 
25 autocorrelation of the source intensity distribution s[k s ) with the vector impulse response 

function h[k\k s ). A generalized bilinear kernel can be formed to include mask blur, as in 
Equation 41 . Resist blur may be included in a generalized bilinear kernel, for example, as shown 
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in Equation 45, where averaging over the depth of the resist stack and/or a range of focus 
positions may be performed. Preferably, a SOCS decomposition of the generalized bilinear 
kernel is performed (Block 540), as discussed further below. The generalized bilinear kernel V 
is formed so as to be independent of the mask transmission function ra(r), which is provided in 
5 Block 550. The image l(r ) at a point r is then computed (Block 560) by combining, or more 

particularly, performing a bilinear integration of the generalized bilinear kernel with the mask 
transmission function, as, for example, in Equation 28, Equation 32, Equation 40, or Equation 
44. The resulting image /(r) can then be used to perform model-based optical proximity 
correction (MBOPC) at Block 570. 

10 In a preferred embodiment of the present invention, the SOCS decomposition of the 

generalized bilinear kernel (Block 540) may be performed as illustrated in FIG. 5. First the 
image plane within the ROI integration domain is gridded (Block 541). If possible, in 
accordance with the invention, the ROI integration domain is preferably folded according to the 
symmetry of the system (Block 542). The generalized bilinear kernel (GBK) is computed and 

15 tabulated at the grid points of the image plane (Block 543), in the folded ROI integration domain. 

In this preferred embodiment, the tabulated GBK values are remapped to a reduced basis (Block 
544), and then the eigenfunctions of the GBK are calculated in the reduced basis, and converted 
back to the original grid (Block 545). If necessary, the dominant eigenfunctions may be 
iteratively refined against the tabulated GBK values, for example, as in the Lanczos method 

20 (Block 546). Then, the convolutions of the dominant eigenfunctions with the possible mask 

polygon sectors are pre-computed (Block 547). 

A preferred embodiment for computing the image (Block 560) is illustrated in FIG. 6. 
First, the mask transmission function that was provided in Block 550 is decomposed into an 
appropriate set of the possible polygon sectors (Block 561). Then, for each of the dominant 

25 eigenfunctions (Block 562), a pre-image is computed by forming the coherent sum of the 

contributions from the pre-computed convolutions of the dominant eigenfunctions with the mask 
polygon sectors within the ROI (Block 563). The pre-image is weighted, preferably by the 
eigenvalue of the eigenfunction, and in particular, more preferably by the square root of the 
eigenvalue. Alternatively, the weights may be determined empirically. Finally, the image l(r) 
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is formed (Block 564) from the incoherent sum of the weighted pre-images of all the dominant 
eigenfunctions. 

The present invention has the advantage that non-scalar effects such as vector effects 
including tailored source polarization, lens birefringence, and resist stack polarization, as well as 
5 blur in the mask or resist, can be incorporated efficiently in image calculation using techniques 

such as SOCS decomposition. For example, the computation of the image (Block 550) can be 
performed efficiently by a method such as outlined in FIG. 5. The mask polygons may be 
decomposed, for example, as sectors (Block 551 ), as indicated by Equation 71 . Similarly, the 
image plane may also be decomposed, for example, by gridding (Block 552). A SOCS 

10 decomposition of the generalized bilinear kernel of the present invention may be performed, 

independently of the mask transmission function (Block 553). The dominant eigenfunctions and 
eigenvalues may be obtained (Block 553), for example, by using the Lanzos method. 

The method of the present invention for computing an image of an integrated circuit 
design may be implemented by a computer program or software incorporating the process steps 

15 and instructions described above in otherwise conventional program code and stored on an 

electronic design automation (EDA) tool or an otherwise conventional program storage device. 
These instructions include providing an impulse response function for the projection lens (which 
may include resist stack effects and blur), the generation of a generalized bilinear kernel, taking 
into account vector effects, and performing a SOCS-like decomposition of the generalized 

20 bilinear kernel, using methods as described above. The instructions to perform the method of the 

present invention may be incorporated into program code to perform model-based OPC. As 
shown in FIG. 7, the program code, as well as any input information required, may be stored in 
EDA tool or computer 440 on program storage device 442, such as a semiconductor chip, a read- 
only memory, magnetic media such as a diskette or computer hard drive, or optical media such 

25 as a CD or DVD ROM. Computer system 440 has a microprocessor 444 for reading and 

executing the stored program code in device 442 in the manner described above. 

It will be appreciated by those skilled in the art that the method and system for 
performing the method in accordance with the present invention is not limited to the 
embodiments discussed above. Accordingly, the invention is intended to encompass all such 
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alternatives, modifications and variations which fall within the scope and spirit of the invention 
and the following claims. 
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